calibration_version = '20191107_nochoice';
file_in = ['../../../data/generated/matlab/simulation/', ...
           calibration_version, '_simulation_fix_capital.csv'];

data = tools.csv2struct(file_in);
const = tools.csv2struct('../../../data/raw/constants.csv');

var.q_B = data.q_B(1);
var.K_B = data.K_B(1);
var.avg_inc_R = data.L_R(1) * data.Y_R(1) / data.emp_R(1);
var.avg_inc_M = data.L_M(1) * data.Y_M(1) / data.emp_M(1);
var.f_M = data.pop_M(1);
var.f_C = data.pop_C(1);
var.w_R = data.avg_wage_R(1);
var.w_C = data.avg_wage_C(1);

psi_M = data.avg_welfare_weight_part_norm_M(1);
psi_C = data.avg_welfare_weight_part_norm_C(1);

var.epsilon = const.epsilon;

% compute reference values for elasticities
wM_wR_0 = data.avg_wage_M(1)/data.avg_wage_R(1);
wM_wR_1 = data.avg_wage_M(end)/data.avg_wage_R(end);
wC_wR_0 = data.avg_wage_C(1)/data.avg_wage_R(1);
wC_wR_1 = data.avg_wage_C(end)/data.avg_wage_R(end);
rel_change_wM_wR = wM_wR_1/wM_wR_0 - 1;
rel_change_wC_wR = wC_wR_1/wC_wR_0 - 1;
rel_change_K_B = data.K_B(end) / data.K_B(1) - 1;

eps_w_M_reference = rel_change_wM_wR/rel_change_K_B;
eps_w_C_reference = rel_change_wC_wR/rel_change_K_B;

tau_B_reference = simple_tax_formula(eps_w_M_reference, eps_w_C_reference, psi_M, psi_C, var)

%% vary elasticities
eps_w_M_range = -5:0.1:5;
eps_w_C_range = -5:0.1:5;
[EPS_W_M, EPS_W_C] = meshgrid(eps_w_M_range, eps_w_C_range);


TAU_B_EPS = simple_tax_formula(EPS_W_M, EPS_W_C,  psi_M, psi_C, var);

figure(1)
contourf(EPS_W_M, EPS_W_C, TAU_B_EPS, 100, 'LineColor', 'none')

% need to transpose TAU_B for proper display in R
TAU_B_EPS_TRANS = TAU_B_EPS';

base_dir = '../../../data/generated/matlab/robustness/';

writematrix(eps_w_M_range', [base_dir, 'EPS_W_M.csv'])
writematrix(eps_w_C_range', [base_dir, 'EPS_W_C.csv'])
writematrix(TAU_B_EPS_TRANS(:), [base_dir, 'TAU_B_EPS.csv'])

%% vary welfare weights
psi_M_range = 1:0.1:10;
psi_C_range = 0:0.1:1;

[PSI_M, PSI_C] = meshgrid(psi_M_range, psi_C_range);


TAU_B_PSI = simple_tax_formula(eps_w_M_reference, eps_w_C_reference, PSI_M, PSI_C, var);

figure(2)
contourf(PSI_M, PSI_C, TAU_B_PSI, 100, 'LineColor', 'none')

% need to transpose TAU_B for proper display in R
TAU_B_PSI_TRANS = TAU_B_PSI';

base_dir = '../../../data/generated/matlab/robustness/';

writematrix(psi_M_range', [base_dir, 'PSI_M.csv'])
writematrix(psi_C_range', [base_dir, 'PSI_C.csv'])
writematrix(TAU_B_PSI_TRANS(:), [base_dir, 'TAU_B_PSI.csv'])